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CNj ■ Abstract 

' A new procedure based on layered feed-forward neural networks for the microplane 

r T ^ material model parameters identification is proposed in the present paper. Novelties 

are usage of the Latin Hypercube Sampling method for the generation of training sets, 
a systematic employment of stochastic sensitivity analysis and a genetic algorithm- 
based training of a neural network by an evolutionary algorithm. Advantages and 
disadvantages of this approach together with possible extensions are thoroughly dis- 
cussed and analyzed. 

q ; 1 Introduction 

Concrete is one of the most frequently used materials in Civil Engineering. Nevertheless, 
as a highly heterogeneous material, it shows very complex non-linear behavior, which is 
extremely difficult to describe by a sound constitutive law. As a consequence, numerical 
simulation of response of complex concrete structures still remains a very challenging and 
demanding topic in engineering computational modeling. 
CnJ ■ One of the most promising approaches to modeling of concrete behavior is based on 

the microplane concept, see, e.g., [71 Chapter 25] for general exposition and [TJ for the 
most recent version of this family of models. It leads a fully three-dimensional mate- 
rial law that incorporates tensional and compressive softening, damage of the material, 
supports different combinations of loading, unloading and cyclic loading along with the 
development of damage-induced anisotropy of the material. As a result, the M4 variant 
of the microplane model introduced in [2] is fully capable of predicting behavior of real- 
world concrete structures once provided with proper input data, see [U |9] for concrete 
engineering examples. The major disadvantages of this model are, however, a large num- 
ber of phenomenological material parameters and a high computational cost associated 
with structural analysis even in a parallel implementation |lUj . Although the authors of 
the model proposed a heuristic calibration procedure [2, Part II], it is based on the trial- 
and-error method and provides only a rude guide for determination of selected material 
parameters. Therefore, a reliable and inexpensive procedure for the identification of these 
parameters is on demand. 

In the view of potential improvements demonstrated in the recent work by Novak 
and Lehky [T2], the applicability of a novel procedure based on artificial neural net- 
works (ANN's) for the microplane parameter identification is examined in the present 
contribution. Individual steps of the identification procedure involve (see also [12} [T3] for 
more details) 
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Step 1 Setup of a, virtual and/or real experimental test used for the identification proce- 
dure. 

Step 2 Formulation of an appropriate computational model. Input data to the model 
coincide with the parameters to be identified. 

Step 3 Randomization of input parameters. Input data are typically assumed to be 
random variables uniformly distributed on a given interval. 

Step 4 Stochastic sensitivity analysis using the Monte Carlo-based simulation. This pro- 
vides us with relevant model parameters which can be reliably identified from the 
computational simulation. 

Step 5 Definition of topology of an ANN used for the identification procedure. 

Step 6 Training of the ANN. The training set is formed from the data generated during 
the sensitivity analysis step. 

Step 7 Validation of the ANN with respect to the computational model. This step is 
usually performed by comparing the prediction of the ANN with an independent set 
of input dataH- 

Step 8 The identification of relevant model parameters using trained ANN from available 
experimental data. 

In the rest of the paper a more detailed description of the individual steps when applied 
to the microplane M4 model is presented. The basic outline of the constitutive model 
together with numerical solution scheme are presented in Section 2. The randomization of 
input parameters using small-sample Monte Carlo simulation is described in Section 3 to- 
gether with stochastic sensitivity analysis. Section 4 is devoted to the procedure of ANN's 
training. The application to the identification of microplane model parameters is intro- 
duced in Section 5. Finally, the results obtained using the methodology are summarized 
in Section 6 together with comments on possible improvements. 

2 Microplane model M4 for concrete 

In contrary to traditional approaches to constitutive modeling, which build on description 
via second-order strain and stress tensors at individual points in the (x, y, z) coordinate 
system, the microplane approach builds the descriptions on planes of arbitrary spatial 
orientations - so-called microplanes, related to a macroscopic point, see Figure [TJ This 
allows to formulate constitutive equations in terms of stress and strain vectors in the 
coordinate system (1, m, n) associated with a microplane oriented by a normal vector 
n. The general procedure of evaluation of a strain-driven microplane model response 
for a given "macroscopic" strain tensor e(x) can be described as follows: (i) for a given 
microplane orientation n normal "macroscopic" strain tensor e(x) is projected onto the 
normal "microstrain" vector e(n) and the shear microstrains e(m) and e(l), (ii) the normal 
and shear microstresses cr(n),<r(m) and <r(l) are evaluated using microplane constitutive 
relations, (iii) the "macroscopic" stress tensor tr(x) is reconstructed from the microscopic 
ones using the principle of virtual work, see, e.g., Chapter 25] for more details. In 
the particular implementation, 28 microplanes with a pre-defined orientation on the unit 
hemisphere is used to evaluate the response of the model. 

1 usually called a "training set" 
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Figure 1: Concept of microplane modeling 



To close the microplane model description, the appropriate microplane constitutive 
relation must be provided to realistically describe material behavior. The model examined 
in the current work is the microplane model M4 [2] . The model uses volumetric-deviatoric 
split of the normal components of the stress and strain vectors, treats independently shear 
components of a microplane and introduces the concept of "boundary curves" to limit 
unrealistically high values predicted by earlier version of the model. As a result, the 
strain-to-stress map e(x) i— > <x(x) is no longer smooth, which complicates the formulation 
of consistent tangent stiffness matrix [TU] and, subsequently, gradient-based approaches to 
material model parameters identification. 

In overall, the microplane model M4 needs eight parameters to describe a certain type 
of concrete, namely: Young's modulus E, Poisson's ratio v, and other six parameters (fci, 

k%, &4, C3, C20), which do not have a simple physical interpretation, and therefore it 
is difficult to determine their values from experiments. The only information available in 
the open literature are the bounds shown in the Table [TJ 
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In the present work, the computational model of a structure is provided by the object- 
oriented C++ finite element code OOFEM 1.5 [13j . Spatial discretization is performed 
using linear brick elements with eight integration points. The arc-length method with 
elastic stiffness matrix is used to determine the load-displacement curve related to the 
analyzed experiment. 



Parameter 



Bounds 



Table 1: Bounds for the microplane model parameters 
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3 Input parameter randomization and stochastic sensitivity 
analysis 



The novelty of the identification approach proposed in [12] is a systematic use of small- 
sample Monte Carlo simulation method for generation of neural network training sets as 
well as stochastic sensitivity analysis. In the particular case of the M4 microplane model, 
each input parameter is assumed to be uniformly distributed on an interval specified in 
Table H 

The Latin Hypercube Sampling (LHS) method [6] is used to generate particular real- 
ization of input variables as it enables to minimize the amount of simulations needed to 
reliably train a neural network. Moreover, the Simulated Annealing optimization method 
available in the software package FREET [11] is used to maximize the statistical indepen- 
dence among individual samples. The Pearson product moment correlation coefficient, 
defined as 

cor = (1) 

vTX^ - x ) 2 YJjji - y) 2 

where x and y denote the expected values of random variables X and Y, is used as 
a sensitivity measure to investigate the influence of individual parameters to a structural 
response. Note that the correlation coefficient is normalized as — 1 < cor < 1, where 
higher absolute values indicate statistical dependence of the random output variable Y on 
the random input variable X. 



4 Artificial neural network 

In this work, layered fully connected feed-forward neural networks with bias neurons (see, 
e.g, [15]) are used for the parameter identification. In general, a neural network is created 
to map the input vector I = (Iq, I\, . . . I m ) on a target vector T = (To, T±, . . . T n ). There 
are L layers denoted as lo,h ■ ■ - Il-i, where Iq is the input layer and Il-x is the output 
layer. The i-th layer Zj has iVj neurons denoted as n-n, n^, ■ ■ ■ n>i,Ni- Each layer except the 
output layer has the bias neuron m t Q. The connections are described by the weights w^ij, 
where Z = 1, 2 ... L — 1 denotes a layer, i = 0, 1 . . . iVj_i is the index number of a neuron 
in the preceding layer Z — 1 (i = for bias neurons) and j = 1,2 ... Ni is the index number 
of a neuron in the layer Z. The output of the neuron nij is then defined as 

Oij = /act f ^0z-i,i.«>J,i,ij , 1 = 1,2. ..L-l, j = 1,2. ..N t , (2) 

O j = Ij, j = l,2...N , (3) 
O lfi = 1, Z = 0,1... L-l, (4) 

where / ac t is an activation function. In our current implementation the activation function 
has the following form: 

■^ act ( S ) = (l + e-a/E) ' ( 5 ) 

where a is the gain of the / ac t- The value a = 0.5 is used in all reported calculations. The 
output vector of each layer Zj is denoted as = (Oj,i, Oj,2 5 • • • O^jyJ. Finally, the neural 
network is propagated as follows: 
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1. Let I = 1. 



2. Calculate M for i = 1, 2 . . . iVj. 



3. / = / + !. 



4. If I < L go to [21 else Ol_i is the network's approximation of T. 



The output error, which is used as a measure of a training level, is defined as 



e = 



m - o L ^) 2 . 



(6) 



i=l 



4.1 Training algorithm 

Behavior of a neural network is determined by a preceding training process. It consists 
of finding the synaptic weights, which have influence on the response of a neural network, 
depending on the different components of an input signal. The training of a neural network 
itself could be considered as an optimization process, because it can be seen as a minimiza- 
tion of neural network output error ([6]) . Then the synaptic weights of a neural network act 
as variables of the optimization algorithm's fitness function. As it was shown earlier, e.g. 
in [3], evolutionary algorithm-based optimizers can significantly outperform the traditional 
methods, e.g. the backpropagation method. Therefore, the evolutionary optimization al- 
gorithm GRADE, proposed in [5], is used for the neural network training. Due to size 
limitations, we present only a sketchy description of the corresponding procedure and refer 
an interested reader to (5J H] for a more elaborate discussion. 

In the computations to follow we will work with the population of 10 x n chromosomes, 
where n is the total number of unknowns in the problem. This population evolves through 
the following operations: 

Mutation Let Xi(g) be the i-th chromosome in a generation g, 



where n is the number of variables of the objective function. If a certain chromosome 
Xi(g) is chosen to be mutated, a random chromosome RP is generated from the 
definition domain and a new one Xk(g + 1) is computed using the following relation: 



Parameter MR is chosen randomly from the interval (0, 1). The number of new 
chromosomes created by the mutation operator is defined by 'radioactivity', which 
is a parameter of the algorithm, with a constant value set to 0.2 for all reported 
calculations. 

Gradient cross-over The aim of the cross-over operator is to create as many new chro- 
mosomes as there were in the last generation. The operator creates new chromosome 
Xi(g+1) according to the following sequential scheme: choose randomly two chromo- 
somes x q (g) and x r (g), compute their difference vector, multiply it by a coefficient 
CR and add it to the better one of them, i.e, 



Xj(fiO = (xa(g),x i2 (g), ...,x in (g)) 



(7) 




(8) 



Xi(g + 1) = max(x q (g);x r (g)) + SG x CR(x q (g) 



Xr(ff)). 



(9) 



5 



The parameter CR is chosen randomly from the interval (0, CL), where CL is a pa- 
rameter of the algorithm equal to 1.0 for all our calculations. SG denotes the sign 
change parameter which is supposed to get the correct orientation of the increase 
with respect to the gradient of the objective function. 

Selection represents the kernel of each genetic algorithm. The goal is to provide a pro- 
gressive improvement of the whole population, which is achieved by reducing the 
number of the "living" chromosomes together with conservation of the better ones. 
Modified tournament strategy is used for this purpose: two chromosomes are cho- 
sen randomly from a population, they are compared and the worse of them is cast 
off. This conserves population diversity thanks to a good chance of survival even 
for badly performing chromosomes. 

Moreover, sometimes the method is "caught" in a local extreme and has no chance to 
escape unless a mutation randomly finds a sub-area with better values. If the gradient 
optimization methods are applied, this case is usually resolved by so-called multi-start 
principle. It consists of restarting the algorithm many times with different starting points. 
Similarly, any type of an evolutionary algorithm could be restarted many times. Neverthe- 
less, the experience shows that there are functions with so-called deceptive behavior (and 
the training of a neural network is one of them), characterized by a high probability that 
the restarted algorithm would fall again into the same local extreme rather than explore 
another sub-area. 

As a solution, the method has been introduced in [4j. It produces areas of 

higher level of "radioactivity" in the neighborhood of all previously found local extremes 
by substantially increasing the mutation probability in these areas (this probability is 
set to 100% hereafter). The diameter of the radioactivity area (finally it defines a n- 
dimensional hyper-ellipsoid for all variables) is set to a 75% percentage of an appropriate 
variable interval. The time of stagnation that precedes the markup of a local extreme 
is a parameter of the method set to 100 generations. Similarly to the living nature, the 
radioactivity in the CERAF method is not constant in time but decreases during the time 
as the solutions produced by the cross-over are trying to get inside the radioactivity area, 
see again [4] for more details. When the number of stagnating generations is determined 
(the change between two best solutions in following generations is less than some very 
small value, 5 • 10~ n in our case), the actual best solution is declared as the center of the 
new radioactive area and the whole population is restarted. 

5 Identification of microplane model M4 parameters 

The present section summarizes the individual steps of M4 material model identification. 
Following the heuristic calibration procedure suggested in [21 Part II], we examine three 
specific experimental tests: (i) uniaxial compression, (ii) hydrostatic test and (iii) triaxial 
test. Advantage of these tests is their simplicity and availability in most experimental 
facilities. Moreover, authors in [2] claim that these experiments are sufficient to determine 
all parameters of the microplane model M4. The results presented in this section can be 
understood verification of this claim. 



Abbreviation of the French expression CEntre RAdioactiF - the radioactivity center. 
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Figure 2: Uniaxial test - Experiment setup 




Figure 3: Bundle of simulated stress-strain curves for uniaxial compression test 
5.1 Uniaxial compression test 

The most common experiment used for the determination of concrete parameters is the 
uniaxial compression test on cylindrical concrete specimens. In particular, the cylinder 
with a radius equal to 75 mm and the height of 300 mm is used. The set-up of the 
experiment is shown in Figure [2 

The LHS sampling procedure has been used to determine the set of 30 simulations re- 
sulting in a "bundle" of stress-strain curves shown in Figure [3j The evolution of stochastic 
sensitivity during the loading process is depicted in Figure 01 The results indicate that 
the most sensitive parameters are Young's modulus E, the coefficient ki, Poisson's ratio v 
(especially for the initial stages of loading) and, for the later stages of loading, the coeffi- 
cient C2o- Therefore, one can expect that only these parameters can be reliably identified 
from this test. 

Moreover, the impact of individual parameters on a position of a peak of stress-strain 
curves is computed. The results of a sensitivity analysis using Pearson's product moment 
correlation coefficient of peak coordinates [e,cx] are listed in Table [21 Results indicate 
particularly strong influence of the k\ parameter, which hopefully allows its reliable de- 
termination. 

Based on the results of sensitivity analysis, the neural network training can be per- 
formed using a nested strategy. First, Young's modulus E with sensitivity ~ 1 in the initial 
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Figure 4: Sensitivity evolution for uniaxial compression test 





Pearson's 


coefficients 


Parameter 


e 


a 


fcl 


0.968 


0.709 


k 2 


0.025 


0.008 




0.015 


0.030 




0.021 


0.074 




-0.019 


-0.020 


C20 


0.158 


0.041 


E 


0.004 


0.684 


V 


0.129 


0.000 



Table 2: Pearson's coefficient as a sensitivity measure of individual parameters to the peak 
coordinates [e,a] of stress-strain curves 

stage is easily identified. To this end, a three-layer ANN is used. In the first layer only the 
neurons corresponding to the values of stresses a z i, cr 2j 2 and a z ^ in the first three points 
of axial strain with extremal Pearson's correlation coefficient are chosen. The second layer 
contains two neurons only; the last layer consists of one neuron corresponding to the pre- 
dicted value of Young's modulus E. For the ANN training the GRADE algorithm is used 
and calculation is stopped after 1,000,000 iterations of the algorithm. 

The three-layer ANN trained using the GRADE algorithm is also used for the identi- 
fication of other microplane model parameters. In Table network's architectures and the 
choice of input values for the identification of each microplane parameter are presented. 



Parameter 


ANN' 


s 


layout 


Input values 


h 


4 - 


2 


- 1 


&x,2, &z, P eak, £z, P eak, prediction of E 




4 - 


2 


- 1 


o" 2 ,20, cr z ,peak, (-z, P eak, prediction of E 


h 


4 - 


2 


- 1 


&z,20, &z,peak, £z, P eak, prediction of E 




4 - 


2 


- 1 


o- z ,20, crz,peak, (-z.peak, prediction of E 




4 - 


2 


- 1 


<7z,20, &z, P eak, ^z, P eak, prediction of E 


C20 


4 - 


3 


- 1 


Cx,iO0; °"z,20) prediction of E, prediction of k\ 


E 


3 - 


2 


- 1 


02,1> C«,2) a z,3 


V 


4 - 


3 


- 1 


&x,i, &x,2, prediction of E, prediction of ki 



Table 3: Neural network architectures 



8 



The results of the identification for an independent set of ten stress-strain curves using 
the proposed strategy are shown in Tabled 



£1 T ' rnpf' pr 
X ClJL dlllctcl 


Absolute error 

average maximal 


Relative error [%] 

average maximal 


h 


2.058e-06 


4.678e-06 


1.34 


2.76 




138.9 


318.7 


38.92 


179.62 


h 


2.679 


5.283 


33.96 


102.33 




52.70 


91.70 


48.33 


107.17 




1.675 


2.278 


37.66 


47.09 


C20 


0.7547 


1.4168 


26.70 


56.69 


E 


229.3 


594.5 


0.74 


1.79 


V 


0.006447 


0.010361 


2.93 


4.72 



Table 4: Errors in the estimated parameters obtained from ten independent tests 

Note that obtained errors are in a close agreement with the results of sensitivity anal- 
ysis. Except E, v and k\, the parameters of the model are identified with very high error 
values. Therefore, additional simulations are needed to obtain these values reliably. 

30 
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Figure 5: Bundle of simulated stress-strain curves for uniaxial compression with fixed val- 
ues of Young's modulus, Poisson's ratio and k\ parameter and one (bold black) measured 
stress-strain curve 

At this point we have to fix already well-identified parameters to the optimized values 
and perform simulations for the rest of parameters. To minimize computational time, 
values for one uniaxial measurement presented later in Section 15.41 were used. Corre- 
sponding values predicted by previously learned neural networks were Young's modulus 
E = 32035.5 MPa and k x = 0.000089046. Poisson's ratio is set to v = 0.2 as a usual 
value of a wide range of concretes. Next, 40 new simulations varying the rest of unknown 
parameters are computed. From this suite, only 34 solutions are valid, i.e. these solutions 
were able to reach the desired value of axial strain e = 0.008. The bundle of resulting 
curves is shown in FigureEJ Note that the black bold curve represents measured data. The 
evolution of Pearson's correlation coefficient during the experiment is shown in Figure El 

The sensitivity shows very high influence of k% parameter at the beginning of the 
loading. If we inspect Figure [71 it is clear that the k 2 parameter influences the stress-strain 
curve only on a very narrow interval and hence the parameter k% cannot be identified from 
this test (the authors of the microplane model indeed proposed k 2 to be estimated from 
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Figure 7: parameter as a function of the stress a 12 (corresponding to e = 0.0011) 



the triaxial compression test). We are more interested in fitting data in the post-peak 
part. For post-peak curves, sensitivity analysis shows especially growing influence of C20 
parameter. This can be demonstrated by a relation between the C20 parameter and a value 
of a stress (crgi) at the end of our simulations, where the correlation coefficient reaches the 
value 0.904429. This relation is graphically illustrated in Figure 

It is clearly visible that this relation is not highly non-linear and any simple regression 
of this data is not possible. We applied the ANN with 3 input neurons chosen to get the 
best prediction of parameters based on post-peak curves. Therefore, one input value is 
a stress value at the peak Upeak 

and the other two inputs are stress values <jqi and ogi 
corresponding to strains e = 0.006 and e = 0.008, respectively. Two neurons in the hidden 
layer were used. Quality of ANN prediction is demonstrated in Figure GO In particular, 
the exact prediction of the searched value corresponds to a point lying on the the axis 
of the first quadrant. Values of predicted parameters are normalized here to the interval 
(0.15,0.85). Dashed parallel lines bound a 5% relative error related to the size of the 
parameter's interval. Clearly, the identification procedure works with an error which does 
not exceed the selected 5% tolerance. 

The attention is also paid to the over-training of the ANN. To control this aspect, the 
evolution of errors in ANN's predictions during the training process on the training and 
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Figure 8: The C20 parameter as a function of a stress (usi) at the end of simulations 




Figure 9: Quality of ANN predictions of C20 parameter 

testing data are monitored (see Figure fTOj) . Recall that if the errors on the testing set are 
much higher than on the training set, we suppose such an ANN to be over-trained. Even 
though this seems to be the case of the current ANN, we attribute such a behavior to the 
fact that there are more training data (in our case 25) then 11 neural weights optimized 
by the algorithm. Also note the typical restarting of the optimization process caused by 
the multi-modal optimization strategy CERAF presented in Section 14.11 

5.2 Hydrostatic test 

The next independent test used for the identification problem is the hydrostatic compres- 
sion test, where a concrete cylinder is subjected to an increasing uniform pressure. Axial 
and two lateral deformations (in mutually perpendicular directions) are measured. The 
experimental setup is shown in Figure ITTb-b. Again, to improve identification precision, 
the parameters E, v and k\ are supposed to be fixed to the previously identified values. 
The "bundle" of stress-strain curves obtained using the LHS sampling for 70 samples is 
depicted in Figure [TTb and the corresponding sensitivity evolution in Figure [T2"l Note that 
the maximal value of a hydrostatic pressure for all these tests is 427.5 MPa. 

The sensitivity information reveal that this test can be used to identify parameter ^3 
from the loading branch while a combination of loading and/or unloading data can be 
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Figure 10: Evolution of ANN's errors durin 






(a) 



(b) 



(c) 



Figure 11: Hydrostatic test, (a) Experiment setup, (b) Initial and deformed finite element 
mesh, (c) Stress-strain curves 



used for foj parameter identification. Moreover, the correlation between the strain at the 
peak of curves and k^ parameter is so high that one can expect their relation to be almost 
linear. This is, however, not the case as illustrated by Figure [131 showing the value of k& 
parameter as a function of a strain s. In spite of a high value of the correlation coefficient 
equal to 0.958586, the noise of these data seems to be very high and a use of a linear 
regression introduces a high error in parameter prediction. 

To identify the k% parameter, which shows high correlation near to the end of elastic 
stage, we evaluate the correlation coefficient between this parameter and a position of the 
end of an elastic stage, k^ and e y i e id correlation is 0.95078 and k% and o y i e id = 0.951873. 
In Figure UM values of k% parameter as a function of these coordinates are presented. 
The discrete values of e y i e id are caused by the size of a time step at the beginning of 
simulations. The noise in data is very high again and it is not possible to reliably use 
a linear regression. Because all other parameters have very small correlation, we suppose 
that the noise of the parameter k-$ is caused by k^ parameter and vice-versa. In other 
words, these noises could be caused by some level of correlation between the parameters 
&3 and k±. Hence we decided to apply an artificial layered neural network again. The 
first 60 simulations prepared by the LHS method were used for training and remaining 
(randomly chosen) 10 simulations for testing. Particular choice of input values as well as 
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Figure 12: Evolution of Pearson's correlation coefficient during the hydrostatic compres- 
sion test for loading (left) and unloading (right) branch 
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Figure 13: k^ parameter as a function of a strain of a peak 



architectures of ANN's is shown in Table To eliminate unknown correlation between 
parameters k^ and k^, their values are used also as inputs into ANN's. 



Parameter 


ANN's layout 


Input values 


^3 


5-2-1 


fell (-yields ^load,2i ^load,5^peak 


ki 


3-2-1 


^3 j Epeafei ^unload,i 



Table 5: Neural network architectures for hydrostatic test 



The architectures of ANN's were chosen manually to get the best precision in predic- 
tions and also to avoid the over-training of the ANN's. Therefore, it is possible to show 
the evolution of ANN's errors (see Figure I15|) during the training. A training process 
with 5000000 iterations takes approximately 20 minutes. Quality of ANN's predictions is 
demonstrated in Figure fTBl Values of predicted parameters are again normalized into the 
interval (0.15,0.85). 

In this way, two ANN's or, in other words, two implicit functions are prescribed. One 
defines a value of k% parameter depending on a value of k^ parameter and some other 
properties of a stress-strain curve, the second define a value of k^ parameter depending 
on a value of k^ parameter and some other properties of a stress-strain curve. Once 
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Figure 14: A:3 parameter as a function of a position of the end of an elastic stage 





Figure 15: Evolution of ANN's errors during the training process in prediction of (a) k^ 
parameter and (b) k^ parameter 



we get some "measured" data and we fix all properties of a stress-strain curve, we get 
a system of two non- linear equations for k^ and k^. We can solve this system, e.g., 
graphically. Both relations are depicted in Figure [T7| for one independent stress-strain 
curve (&3 = 7.84293, k& = 155.551). Their intersection defines searched parameters, 
A)3 = 8.15687 and fej = 154.072 in this particular case. The precision of the proposed 
strategy is visible in comparison of corresponding stress-strain curves, see Figure [TH1 Note, 
that E, v and k% were the same as in previous section and remaining parameters, i.e. &2, 
C3 and C20, were chosen randomly]! 



5.3 Triaxial test 

The last experiment, used for the purpose of parameter identification, is a triaxial com- 
pression test. To this end, a specimen is subjected to the hydrostatic pressure an- After 
the peak value of an is reached, the axial stress is proportionally increased. The "excess" 
axial strain e = et — £h, where et and ejj denote the total and hydrostatic axial strain, 
is measured as a function of the overall stress a. The test setup is shown in Figure [T9l 

3 Theoretically, the value of C20 is known in this stage from the previous identification step. Nevertheless, 
the numerical simulation of the unidirectional experiment takes substantially more time to complete, see 
Section [H] for an example, and therefore the independence of C20 allows us to proceed with the inverse 
analysis even though the first phase is not finished. 
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(a) (b) 
Figure 16: Quality of ANN prediction of (a) ks parameter and (b) £4 parameter 
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Figure 17: Relations of k% and k^ parameters 

At this point, we assume that parameters E, v, k±, k^ and k^ are known from previous 
identifications^. Next, 70 simulations (60 training and 10 testing) of the triaxial com- 
pression test are computed by varying three remaining parameters k2, C3 and C20 The 
bundle of stress-strain curves for an = 34.5 MPa is shown in Figure [2"U1 together with the 
evolution of Pearson's correlation coefficient during the experiment in Figure I2T1 





Pearson's 


coefficient 


Parameter 


e 


a 


k 2 


0.585 


0.791 


C3 


-0.067 


-0.088 


C20 


0.664 


0.329 



Table 6: Pearson's coefficient as a sensitivity measure of individual parameters to the peak 
coordinates [e,a] of stress-strain curves 

In addition, the correlation coefficient between microplane parameters and stress and 
strain values of peaks is computed. These results are shown in Table El It is visible that 
maximal correlation is between k 2 parameter and the value of the stress 029 corresponding 
to the strain equal to €29 = 0.01276. This correlation is 0.88956 and at the same time, 
the correlation between these 029 values and C20 parameter is very small, therefore C20 
parameter does not influence the relation between &2 parameter and 029- Figure [2"2"1 shows 

4 i.e. E = 32035.5 MPa, v = 0.2, fci = 0.000089046, k 3 = 8.15687 and k 4 = 154.072. 




k4 parameter as a function of k3 parameter 
k3 parameter as a function of k4 parameter 
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Figure 18: Comparison of original simulation and simulation for predicted k% and k± 
parameters 




(a) (b) (c) 

Figure 19: Triaxial compression test, (a) Experiment setup, (b) Initial and deformed mesh 
at the end of hydrostatic loading, (c) Initial and deformed mesh at the end of total loading 



that only small values of C20 parameter disturb this relation. In particular, points related 
to the C20 parameter smaller then 1 lead to oscillatory dependence. 

Because the highest correlation for hi parameter is again at the beginning of the loading 
and after our experiences in identification parameters from the uniaxial compression test, 
we were again afraid of small significance of k<i parameter to the shape of curves. Also 
CT29 does not seem to be significant. Therefore we made several short computations with 
randomly chosen fixed value of C20 parameter to filter out its influence. We have got 
a bundle of curves showing similar spread of values as curves in Figure [2Ql Therefore, it 
can be concluded that these differences are probably caused by k<i parameter only and a 
neural network for hi parameter identification can be designed. Because the bundle of 
curves varies mostly in the post-peak part and we would like to get a predictor capable 
to fit this part of a curve properly, we use a pea k and <7ioo as input values. The latter 
one correspond to the end of our simulations, where e = 0.044. We also add the third 
input value - 029 - because of its small correlation with C20 parameter. Two neurons in 
the hidden layer were used. Quality of the ANN prediction is demonstrated in Figure [231 
Under- and over-fitting issues were again checked by errors evaluations during the training 
process, see Figure [Ml 
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Figure 20: Bundle of simulated stress-strain curves for triaxial compression test 
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Figure 21: Evolution of Pearson's correlation coefficient during the triaxial compression 
test 

Almost perfect precision of the proposed strategy is visible in comparison of corre- 
sponding stress-strain curves for &2 = 748.857 and its prediction equal to 767.777 and 
randomly chosen parameters C20 and C3 , see Figure [25j 

5.4 Application to measured data 

In previous sections, we have shown that the proposed methodology is able to identify all 
but one (03) parameters from computer-simulated curves. To demonstrate the applicability 
of the proposed procedure, a real simulation should be examined. However, only limited 
experimental data from uniaxial compression tests are available to authors which leaves 
us with only one uniaxial stress-strain curve to be identified. As was mentioned previously 
in Section 15.11 Young's modulus E = 32035.5 MPa, Poisson's ratio v = 0.2 and k\ = 
0.000089046 are predicted by the neural network for this measurement. Next, 30 samples 
for parameter C20 are computed, see Figure [5j If we zoom into the loading part of a 
stress-strain curve, it is clear, that the real measurement is under all simulated data, see 
Figure [26l This part is influenced by high correlation of k^ parameter and therefore, 
it is clear that the &2 parameter cannot be obtained from this test. Finally we applied 
our trained ANN to predict the C20 parameter for the measured data and we have got 
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Figure 23: Quality of ANN prediction of hi parameter. 



C20 = 5.27065. This value is out of the interval specified for this parameter, but it is 
not surprising since it is visible in Figure [5] that measured data somewhat deviate from 
the simulated bundle of curves. The final comparison of measured data and a simulation 
for predicted values of E, u, k\ and C20 parameters is shown in Figure [271 The rest of 
unknown parameters are same as in previous sections. 



6 Conclusion and future work 

In the present contribution, an example of the engineering problem, which is difficult to be 
solved by traditional procedures, was solved using methods of soft computing. Particularly, 
an artificial neural network was used to estimate required microplane material model 
parameters. As the training procedure, the genetic algorithm-based method GRADE 
was used. A number of needed simulations is reduced by the application of the Latin 
Hypercube Sampling method. The sensitivity analysis shows not only the influence of 
individual parameters but also approximately predicts the errors produced by the neural 
network. 

Results of our contribution, see Table [7J confirm the claims made by authors [2 J of the 
microplane M4 model on individual parameters fitting. Only the parameter C3 remains 
undetermined but the parameter C3 should be almost constant for the wide range of con- 
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Figure 24: Evolution of ANN's errors during the training in prediction of &2 parameter 
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Figure 25: Comparison of original simulation and simulation for predicted parameters of 
triaxial compression test 

cretes and our computations confirm almost zero impact of this parameter on stress-strain 
curves. 

The rather severe disadvantage of the microplane model, and also of the proposed 
methodology, is an extreme demand of computational time. A suite of 30 uniaxial tests 
consumes approximately 25 days on a single processor PC with the Pentium IV 3400 MHz 
processor and 3 GB RAM. If we run tests in parallel on 7 computers, the needed time 
is less than 4 days. The hydrostatic and triaxial tests are less demanding, by running in 
parallel on 7 computers the required time is less than one day for each test. Although 
several neural networks were created to avoid the time-consuming numerical analysis, 
the proposed methodology still needs to compute 30 uniaxial tests to properly identify C20 
parameter and a set of 30 hydrostatic and triaxial tests to fit k^, k^ and &2- This drawback 
will be the subject of future work. 
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Figure 26: Bundle of simulated stress-strain curves for uniaxial compression and one (bold 
black) measured stress-strain curve under zoom 
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Figure 27: Comparison of measured data and results of final simulation. 
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Parameter 


Test 


ANN's topology 


TP 
Fj 


Uniaxial compression 


9 1 1 1 
+ z + 1 


V 


Uniaxial compression 


4 + 3 + 1 


fcl 


Uniaxial compression 


A 1 O 1 1 

4 + 2 + 1 




Triaxial loading 


3 + 2 + 1 




Hydrostatic loading 


5 + 3 + 1 




Hydrostatic loading 


3 + 2 + 1 




X 


X 


C20 


Uniaxial compression 


3 + 2 + 1 



Table 7: Final status of M4 identification project 
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